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ABSTRACT 



The specific impetus for this work was a conceptual design of a 
submarine using the toroid as the pressure hull. This work is a 
continuation of the work started by Bowen (1987). As such, it is hoped 
that a better understanding of the behavior of a toroid under 
hydrostatic pressure can be realized. 

This work began with a review of efforts to solve complete toroidal 
structures. A specific toroid was then modeled in the B0S0R4 computer 
program to obtain displacements of the meridian under hydrostatic 
pressure. Functions were then derived that described the general form 
of these displacements. Using these functions as the assumed solution 
for the energy method an energy balance was made and a program was 
written to solve for the displacements of a generic toroid under 
hydrostatic loading . 
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CHAPTER 1 
INTRODUCTION 

The purpose of this work was to examine the response of a thin 
circular toroidal shell under external hydrostatic pressure. The 
approach used was to model a specific toroid in the B0S0R4 computer 
program (Bushnell, 1977) and then use the output results for 
displacements to generate functions that modeled these displacements. 
Using the functions obtained that model the displacements of the toroid 
under hydrostatic pressure, an analytical analysis was performed with 
these functions used as the assumed solution for the energy method. The 
energy method generated the equations necessary to solve for the 
displacements of a generic circular toroid. A program was then written 
to solve the equations generated in the energy method so that 
displacements for any specific toroid under hydrostatic pressure could 
be obtained. 

The analysis conducted only considers displacements in the linear 
elastic range. For the shell being analyzed the following assumptions 
were made: 

1. The material of the shell is isotropic and homogeneous. 

2. The thickness of the shell is constant. 

3. The thickness of the shell is small compared to the radii of 

curvature . 

4. The displacements are symmetric about the X-Y plane (see 

Figure 1.1). 

5. Thin shell theory holds: normals to the undeformed surface 
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remain normal. 



For a symmetrically loaded shell, with the above assumption of 
symmetric deformation, a small displacement of any point can be resolved 
into two components: v in the direction of the tangent to the meridian 
and u in the direction of the normal to the middle surface (Timoshenko 
and Woinowsky-Krieger , 1959). Therefore, the only displacements 
discussed throughout the text will be u and v. 

The toroid's geometry is straight forward but in order to be 
consistent throughout the text the following definitions will be made 
and are illustrated in Figure 1.1: 



DEFINITIONS: 

R : Major radius of toroid. Radius of rotation about the Z axis 

of the circle of radius r to form the toroid. 

r : Minor radius of toroid. Radius of the circle which is rotated 

about the Z axis. 

h : Thickness of the shell (assumed to be constant) . 

S : Angle of rotation measured counter clockwise from the X-Y 

plane at a distance R + r from the origin. 

<f> : Angle of rotation about the Z axis measured counter clockwise 

from the positive X axis. 
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Figure 1.1 



Description of Geometry 
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CHAPTER 2 
COMPUTER MODELING 

A computer solution for the displacements of the toroid under 
external hydrostatic pressure was sought to determine the shape of these 
displacements for a specific toroid. The computer program used to model 
the toroid was B0S0R4 (Bushnell, 1977). The B0S0R4 computer program is 
a hybrid finite difference - finite element program used to analyze 
complex shells of revolution. This program was developed by David 
Bushnell at Lockheed Missiles & Space Co., Palo Alto, California. 

Although specific reference was not found that indicated that 
B0S0R4 could adequately model a toroid, it was believed that since 
B0S0R5, a similar program that considers elastic -plastic material 
behavior, had been used to model a toroid that B0S0R4 would also 
adequately model this structure (Bushnell, 1985). 

The specific toroid that was analyzed with B0S0R4 was the toroid 
specified as the default values by Bowen (1987). The toroid geometry 
used was as follows: 

R - 8 inches 
r - 2 inches 
h - 0.02 inches 

Young's Modulus, E - 30 x 10 6 psi 
Poisson's Ratio, v - 0.3 

Modeling the toroid in B0S0R4 was accomplished by taking advantage 
of the symmetry of loading of the structure about the X-Y plane. 
Therefore, only half of the structure was discretized in the program. 
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Because only half of the structure was modeled in the program, boundary 
conditions had to be established at the symmetry points. The boundary 
conditions imposed on the structure are as follows: 



@ 



@ 



6 - 0 radians 

Tangential Displacement, (v) - 
Rotation in the 6 direction, 

6 - II radians 

Tangential Displacement, (v) - 
Rotation in the 6 direction, 





0.0 

- 0.0 



0.0 

- 0.0 



The normal displacement (u) is defined as positive in the direction 
of the outward pointing normal. The tangential displacement is defined 
as positive in the direction of increasing angle 6 (See Figure 2.1). 

Having the toroid modeled in B0S0R4, the program was run to 
determine the pressure required to buckle the structure. The thought 
being that if the buckling pressure predicted by B0S0R4 was consistent 
with the buckling pressures predicted by theory then the pre -buckling 
displacements generated by B0S0R4 could be used to predict the general 
form of the displacements of a generic toroid. 

The buckling pressure obtained from B0S0R4 for the toroid modeled 
differed by less than 10% from the predicted buckling pressure as 
presented by Sobel and Flugge (1967) . Since the buckling pressure 
obtained was consistent with the predicted pressure and the predicted 
buckling pressure was substantiated experimentally by Almroth, Sobel and 
Hunter (1969), it was concluded that the results from B0S0R4 for the 
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Figure 2.1 
B0S0R4 Modeling 
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pre-buckling displacements could be used to predict the general form of 
the displacements of a generic toroid loaded under uniform hydrostatic 
pressure . 

In any finite difference or finite element program an investigation 
into the number and spacing of the nodes used to describe the structure 
must be undertaken to ensure that the structure is adequately described. 

For this investigation, equal spacing of the nodes was used throughout. 

To investigate the effect of node distribution, the number of nodes 
were doubled and thus the spacing was reduced by a factor of two. The 
results for both the magnitude of the buckling pressure and the shape 
and magnitude of the displacements did not change with the increase in 
the number of nodes. It was therefore concluded that the structure was 
adequately described and that the results obtained were valid. 

Figure 2.2 is a graphical presentation of the normal displacement 
(u) and figure 2.3 is a graphical presentation of the tangential 
displacement (v) obtained from the B0S0R4 output for the toroid modeled. 
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Figure 2.3 
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CHAPTER 3 



CURVE FITTING 

The next step in this study was to generate a function that would 
analytically describe the normal and tangential displacements of the 
toroid under hydrostatic loading. 

First a check was made, using the B0S0R4 computer program, to 
ensure the displacements were linearly dependent on the applied 
pressure, as should be the case since this study was only considering 
behavior of the toroid in the linear elastic region. As can be seen 
from figure 3.1 the normal displacement is indeed linearly dependent on 
the applied pressure. Although not included here, the tangential 
displacement is also linearly dependent on the applied pressure. 

To generate a function that would adequately describe the 
displacements, a form of the displacements must be assumed and then a 
curve fitting technique must be used to fit the data to this assumed 
form. 

Since the displacements are considered to be symmetric about the 
X-Y plane and considering the way in which the positive direction for 
the normal displacement (u) and the tangential displacement (v) were 
defined in Chapter 2, the following statements can be made: 

For the normal displacement (u) , 

u(0) - u(-0) (3.1) 

which is the definition of an odd function. 

For the tangential displacement (v) , 

v ( 0 ) - - v ( - 0 ) (3.2) 
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Figure 3.1 
Incremental Loading 



Xormal Displacement 




which is the definition of an even function. 

The functions assumed to describe the displacements were Fourier 
Series for both the normal (u) and tangential (v) displacements. Since 
both the normal and tangential displacements have a period of 2nr (the 
arc length of the circular cross section) and noting the above 
statements about even and odd functions, the Fourier Series for the 
displacements can be written as follows: 

u(0) - + X COS (3.3), 

v(0) - £ bn sin fcp] (3.4). 

n-1 ^ J 

For the displacements to be symmetric about the X-Y plane, or a 
half period of the function, x and L can be defined as follows: 

x - rd , 

L - *r 

From these definitions, u and v can be rewritten as: 

00 

u(0) - + £ an cos (n0) (3.5), 

n-1 

oo 

v(0) “ £ b" sin (nff) (3.6). 

n-1 

While the Fourier series described above will give an exact 
representation of the displacements, the goal of a curve fitting routine 
should be to determine how many terms of the series must be used to 
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adequately describe the function. To determine the number of terms 
required to describe the displacements a least squares curve fitting 
program was written to be used in conjunction with the output data from 
B0S0R4. The programs were adapted from a curve fitting routine 
presented by James, Smith and Wolford (1977) and are included in 
Appendix A. 

Using these programs to determine the coefficients of the Fourier 
Series it was found that five terms were inadequate to describe the 
displacements, as can be seen graphically in Figure 3.2, and it was 
determined that ten terms were required to adequately describe the 
normal displacements, (See Figure 3.3), and nine terms were required to 
describe the tangential displacements. Although only the normal 

displacement (u) is presented in Figures 3.2 and 3.3 similar results 
were obtained for the tangential displacement (v) . 

With the number of terms determined that are required in the 

Fourier series u and v can be written in final form as: 

9 

u(0) - + £ an cos (nff) (3.7), 

n * 1 
9 

v(ff) - l bn sin (n6) (3.8) . 

n ■ 1 

These definitions of u and v will be used in the next section as the 
assumed forms of the displacements of the toroid. 
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NORMAL DISPLACEMENT 



Figure 3.3 
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CHAPTER 4 



ENERGY METHOD 

4 . 1 Introduction : 

To determine the displacements of a structure when it is subjected 
to an external load, the external load must be balanced by the internal 
forces of the structure caused by the external load. 

For this study, the method of minimizing the total potential 
energy, as described by Shames and Dym (1985), was used. The energy 
method was chosen in hopes of avoiding the singularities in the solution 
at 0 - ± - ,(Reissner, 1963 and Senjanovic , 1972). The total potential 
energy is defined as: 

H - U - W (4.1.1) 

where U and W are defined as; 

U : internal energy of the structure, 

W : potential energy of the applied loads. 

To determine the minimum potential energy, and thereby determine the 
displacements caused by the external load, the first variation of the 
total potential energy is set to zero. 

6 (1) n - 0 (4.1.2) 

The first variation of the total potential energy can be defined as the 
following partial differential equation: 




i 



where x, is any general parameter used to describe U and W. 

Since the first variation of the total potential energy is set to zero, 
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the equation to be solved to determine the displacements of the 
structure can be defined as; 



aw _ au 
3x ax 



(4.1.4). 



i i 

In order to solve the above equation for the displacements caused by the 
applied load, U and W must be defined as functions of the displacements . 

As described by Bushnell (1984), the strain energy of the structure 
can be defined as follows: 



U - U + U 
s c 



(4.1.5) 



where ; 

strain energy due to elongation and changes in 
curvature , 

: strain energy due to constraints or boundary conditions. 
For the toroid, the symmetry about the X-Y Plane will be taken 
advantage of and only half of the structure will be analyzed. The 
following additional geometric definitions are required for the analysis 
(See Figure 4.1). 



DEFINITIONS: 

a : Ratio of the major and the minor radius of the toroid. a = - 

r : Distance from the Z axis to the meridian of the toroid. 

1 

r^ - R + r cos(6) - r[ a + cos (0)] 

r : Radius of curvature in the 6 direction. 

2 

r t _ R + r cos(fl) 

2 cos(0) cos(0) 
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Figure 4.1 
Additional Geometry 
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4.2 Strain Energy: 

The strain energy due to elongation and changes in curvature can be 
defined as: 

U - - f a € dv , i.j - e , 4> (4.2.1) . 

S 2 J v ij ij 

The strain energy due to constraints, U , will be addressed later. 

Because this study is only considering hydrostatic loading of the 
toroid, the loading will always be normal to the surface of the toroid. 
Therefore, there will not be any shear stresses or shear resultants from 
the loading imposed. 

a - e - 0 , if i * j . 
ij ij 

The 9 and 4> directions are thus principle directions and the strains can 
be represented by reference surface strains and changes in curvature. 
The strain energy can be defined in terms of forces and moments which 
are defined as an integral of the stresses through the thickness of the 
shell. The strain energy can then be represented by a surface integral 
as follows: 



U - - f f N c + M K 1 ds (4.2.2) . 

S 2 J 5 ^ ij ij ij ij J 



For the case of hydrostatic loading the strain energy becomes: 



U - - f [ fN . £ . + N.£. ) + + M ,K , ] 

s 2js(_'-00 <t> <t> J ^ e e <t> 4 > J 



ds 



(4.2.3) 



9 d 4> 

where the forces and moments are defined as follows (Timoshenko and 



Woinowsky - Krieger, 1959): 

N * “ U 



> [ 1 • i ) dz 

\ • f. % ( 1 • l ] dz 



(4.2.4) , 

(4.2.5) , 
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M « - It °» z c 1 - ! ) dz 



(4.2.6), 



- J. % * [ 1 - f ,] 



dz 



(4.2.7). 



Since the assumption has been made that normals to the undeformed 
middle surface remain normal, the strains as a function of thickness can 
be expressed in terms of reference surface strains and changes in 
curvature (Timoshenko and Woinowsky-Krieger , 1959) . 

- z (4.2.8) 



e 4> t “ U ' Z K 4> 



(4.2.9) 



The terms, e n and e. , represent the total strain in the 6 and 4> 
u <p 

t t 

directions and k a and k represent the changes in curvature of the 
v <p 

reference surface . 

From Hook's Law the stresses can be defined as follows: 

E 



4: ( u + v, i 

•u ) v t t J 

e r 

" 2 . [ U + ve 6 

(l-u)'" t t J 



(4.2.10) 



(4.2.11) 



substituting the expressions for e. and e , , 

V <p 

t t 



O n — 



°<t> " 



(l 



(l 



■ I e A + ll€ , -zf + UK. 1 

, u 2 ) I e <t> e 4> J j 

I £, + ve. -zf K. + l/K„ 1 

. v 2 ) L * 8 L 4 8 J J 



(4.2.12), 



(4.2.13) . 



Taking the middle surface as the reference surface and neglecting 



z z 

the terms, — and — , as small as compared to unity, integration over 
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the thickness of the shell can be represented by: 



„+h/2 



N # - I 

- h 7 2 ( 1 



— f 

-v 2 ) [ 



e e + V€ 4 > ‘ z ( K e + VK * 



3 ) 



dz 



N„ - 



-it- - f . ♦ 1 

<l-v 2 ) l * *> 



+h/2 



N, - J 

* -h J /2 (1 



+ - z ( % + dz 

) 






Eh J" 
L-u 2 ) L 

f + h/2 _ r 

M «-J r ['» + V U <“ 

-h/2 (l-u ) 



+ 



J) 



s + II z dz 



M, 



• Eh r «, + i 

12 (1* u 2 ) J 

r +h/2 e f 1 

* Ji (l-v 2 ) l* ' 1 * 1 J J 



M 



Eh" 



* 12(l-u 2 ) 1 * * 



4- ] I z dz 



( % + •'** ) 



(4.2.14) , 

(4.2.15) , 

(4.2.16) , 

(4.2.17) , 

(4.2.18) , 

(4.2.19) , 

(4.2.20) , 

(4.2.21) . 



Combining the expressions for the forces and the moments and 
substituting these back into the expression for the strain energy, the 
strain energy becomes: 



U 




Eh 

(l-u 2 ) 




+ 2 u v* + 




Eh" 
12(l-u") 



f 2 2 

I k 9 + 2 vk 6 K * + K <(> 



ds 



(4.2.22) . 



4.2.1 Strains: 

The next step is to define the surface strains in terms of the 
displacements u and v. This will be accomplished by following the 
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presentation by Timoshenko and Woinowsky-Krieger , ( 1959 ) . 

Considering an element AB of the meridian, (See Figure 4.2), the 

strain in the d direction can be defined as follows: 

€ due to v , 

9 



€ - 



3v 

v + 39 dd ' V 



r dd 



1 3v 
r dd 



(4.2. 1.1) 



« due to u can be represented by the average change in length of 
9 

the element divided by the original length, 



r dd 



(4. 2.1. 2) . 



The component, - ^ dd , will be neglected as a small quantity of higher 
2 6 9 

order. The total strain in the d direction can then be represented as: 



If ^ 8v } 

-r [ U + 5?J 



(4. 2. 1.3). 



The strain in the 4> direction can be defined as the change in r^ 
due to u and v divided by the original length of the element (See Figure 



4.3). 



A r T - [ u cos(0) - v sin(0) ) d^> , and 

( u cos(0) - v sin(0) ] & 4 > 

*<t> r T d 4> 

« , - — — t . .. . [u cos(^) - v sin(0) 1 

4> r (a + cos (d)) ^ J 



(4. 2. 1.4), 
(4. 2. 1.5). 



4.2.2 Change in Curvature : 

The change in curvature of the middle surface will now be defined. 
Again considering the element AB of the meridian, (See Figure 4.2), 
the rotation in the d direction can be defined as follows: 
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Figure 4.2 

Differencial Element 
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Figure 4.3 



Radial Displacement 



I 
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£ due to v, 



0 




de 



(4.2.2. 1) 



0 due to u, 

— de 

0 — - de (4. 2. 2. 2). 

^ r r 50 

The total rotation in the 6 direction is; 

*, -H V • &] dS ('•.2.2.3, . 

The change in curvature in the 6 direction can be represented by the 
change in rotation in the 6 direction divided by the undeformed arc 
length in the 6 direction. 



1 d ^9 If 5v a 2 u ] 

' r 09 r 2 l « ' a * 2 J 



(4. 2. 2. 4) 



Because the assumption has been made that the deformations will be 
symmetric, the rotation in the 4> direction will be opposite to the 
rotation in the 6 direction in magnitude and rotated into the <f> surface 
by sin (0) . 



0* “ * fi e sin (6) - i( - v ] sin (6) d 4 (4. 2. 2. 5) 

Although the rotation in the 4> direction is constant in <j> , it does vary 
in the 6 direction. Therefore, the change in curvature due to loading 
can be represented by the rotation in the 4> direction divided by the 
undeforraed arc length in the <f> direction. 



r (a + cos(0)) d<f> 



(4. 2. 2. 6) 
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and substituting the expression for 



4 > ' 



sin (0) 



r (a + cos 



( 8 )) f 



— 

89 



- V 



(4. 2. 2. 7) . 



4.2.3 Constraint Energy: 

A complete toroid has no constraint or boundary conditions when 
loaded only by hydrostatic pressure. However, since only half of the 
toroid is being analyzed here, boundary conditions must be imposed at 
the symmetry points, 9 - 0 and 9 - i r. The boundary conditions that are 
required at these two points are: 

The vertical displacement is zero, v - 0, 

and the rotation in the 9 direction is zero, ^ - 0. 

o u 

As Bushnell (1984) points out the constraint energy can be accounted for 
with the introduction of Lagrange multipliers and can be represented as 
follows : 



U c - * \» V <°> + \J ?<”> + V V <'> 



(4.2.3. 1) . 



In this formulation of the constraint energy the Lagrange multipliers 
are unknowns and are solved for along with the displacements when the 
total potential energy is minimized. The actual form of the constraint 
energy will be discussed later. 



4.3 Potential Energy of Applied Loads: 

The potential energy of the applied loads is simply the work done 
by these loads, which can be represented by the applied load multiplied 
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by the deflection in the direction of the applied load. For the case of 
hydrostatic loading, the applied load is always normal to the surface 
and the displacement in the direction of this load is -u for this study. 

Therefore, the potential energy of the applied hydrostatic load can be 
represented as: 

W - - P u ds (4.3.1) 

where P is the external hydrostatic pressure applied and is 
positive in the direction of -u. 



4.4 Minimization of the Total Potential Energy: 

From section 4.1 it was shown that to minimize the total potential 
energy the following equation needed to be solved; 



an 

ax 



o 



(4.4.1), 



i 

and it was shown that the above equation could be represented by 
the following equation: 



aw _ au 
ax Sx 



(4.4.2). 



i i 

To present the form of the above equations the following simplifications 
in notation will be used: 



u' 



au 

39 



u' 




V 1 



3v 

39 



C 



Eh 

(iV) 



Eh“ 



12 ( 1-0 
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4.4.1 Internal Energy: 

Since it has been shown that none of the terms in the internal 
energy depend on <f> the surface integral can be rewritten and integrated 
with respect to 4 > . 






e: + 2 u V* + V ' 



Eh' 



12 (1-u 2 ) * 



+ 2vKgK^ + *Pjr( a + cos(*))d* rd« <4. 4. 1.1) 



Which when integrated with respect to <f> becomes, 



U - jt r 

s 



12 



f [ [ + 2u£.£, + £ 2 ) - 

l l(l-v 2 ) 9 9 + +* 

^"7 + 2vK e K 4, + *p]( a + cos (^)]d« 



(4. U. 1.2) . 



Using the above simplifications in notation and substituting the 
expressions presented in section 4.2.1 and 4.2.2 the strain energy 
becomes : 



U - * r 2 £ ^2 |^ [ u 2 + 2uv' + v' 2 ] + 

(~ a ~ T cos(g)) C (u 2 +uv')cos(^) - (uv + w')sin(tf)) 
+ "os ( ff ) ) 2 C u 2 cos 2 (0) - 2uv cos(0)sin(0)+ 

in 2 (8) ) j - £ [v’ 2 - 2v'u''+u'' 2 ) + 



2 

v s 



2 u sin J c , < , , , , . , 

7 ; TT\\ I V'U - W' -u'u'' +VU' ] 

(a + co s(8 ) ) J 
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(4.4. 1.3) . 



sin ( 0 ) 
(a + co s(9) 



^ 2 [u' 2 - 2u'v + v 2 ] j j (a+cos(0)] dd 



The above equation represents the strain energy, U g . However, for 
minimization of the total potential energy the expression that is needed 

au 



IS 



dx 

i 



£ - * fl C [ C * 2v-g’) 



2v C/o , 3v' , . . 3v 3u , 3v' ,3v . ■> 

7 7T r T [(2 u -T- +U-T— +V't- COS(« • U-T-+V— +V-T — +v'-r- )sin(0)J 

(a + cosCff))'- 3x 3x 3x 3x 3x 3x 3x J 



(a + L«) ) 2 Oh cos2 < # >- 2<u H +v H ><=°s<<Dsin(#> + 2vg S in 2 («)] 

i i i i J 

‘ (*■£’- 2 < v 'i' ; * 



D 

- -2 
r 



2 t/ sin(0) r f du f ,3v' 3v' 3v ,3u'' , ,3u' 3u' ' , , 3v 

(a + cos(0))V 3x 3x 3x 3x 3x 3x 3x 3x J 



(a 



i i i i 



The above equation represents the first term on the right hand side of 
the following equation that will be needed to minimize the total 



potential energy and thereby solve for the displacements, 



aw 

ax 



au 



au 



3x dx 
i i 



(4.4. 1.5) . 



4.4.2 External Energy: 

Similar to the internal energy, the external energy does not depend 
on 4> and the surface integral can be rewritten and integrated with 
respect to <£. 



n 2 7T r * 

P u r[ a + cos (0) ] d<j> r d6 



(4. 4. 2.1) 
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Which when integrated with respect to <f> becomes, 

r* 

W - -2P*r 2 J u ( a + cos(0) ) d 6 (4. 4. 2. 2). 

For minimization of the total potential energy the expression that is 

. . . aw 

needed is — 
dx 

i 

- -2P*r 2 J C a + cos(*) ] d* (4. 4. 2. 3) 

i 0 i 

This completes the expressions needed for the minimization of the total 
potential energy, with the exception of the constraint energy, U 

c 

4.5 Assumed Functions: 

In Chapter 3 it was shown that the displacements u and v could be 
represented by a Fourier series. These functions of 6 will be used as 
the assumed form of the displacements and the coefficients of the 

Fourier series will be solved for from the equations presented in the 
minimization of the total potential energy. The assumed functions for u 
and v are different from those presented in Chapter 3 in that the 
constant term in the expression for u has been included inside the 

summation. 

9 

u(0) ~ X 311 cos ( n ^) (4.5.1), 

n - 0 
9 

v(0) - £bnsin(n0) (4.5.2). 

n -0 

The above expressions for u and v involve 19 unknown coefficients that 
will have to be solved for to describe the displacements . 

Using these expressions for u and v the constraint energy, U 

c 
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discussed in section 4.2.3 can now be addressed. The constraints needed 
were the boundary conditions at the symmetry points, 6 - 0 and 6 - n. 
These boundary conditions are: 

v - 0. Since sin(n0) - 0, for 6 - 0 and sin(n*) - 0, for all n, 
these boundary conditions are satisfied by the assumed form of v 
and additional constraints are not needed. 

g 

4^-0. Since - - Y n an sin (n0) , then — ■ 0 at ^ - 0 and 

Ou 0 u Ou 

n *0 

6 - n for the same reasons as stated above for v. 

Since no additional constraints are needed to meet the boundary 
conditions, Lagrange multipliers are not needed and the constraint 
energy, U - 0. 

4.6 Solution: 

With the assumed forms of the displacements defined, the problem is 
now restricted to solving for the 19 unknown coefficients of the Fourier 
series. Using the following substitutions the equations can be set up 
to solve for these coefficients: 



u - l an cos (n0) 

n ■ 0 

9 


(4.6.1) 


u' - - £ n an sin (n0 ) 

n *0 
9 


(4.6.2) 


u' ' - - Y n a n cos (n0) 

n « 0 
9 


(4.6.3) 


v - Y bn sin (n0) 

n » 0 
9 


(4.6.4) 


v' - £ n bn cos (n0 ) 


(4.6.5) 



n *0 
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Substituting the above expressions into the following equation will then 
represent the 19 equations needed to solve for the 19 unknown 
coefficients : 



aw 

ax 



au 

s 

ax 



(4.6.6) . 



(4.6.7) 



i i 

Which can be represented in matrix notation as follows; 

[■>]-[ F ] [\] 

where the above terms are defined as: 

aw 

D : a 19 x 1 column consisting of the terms, — , 

ox 

i 

au 

F : a 19 x 19 matrix consisting of the terms, — * , 

OX 

i 

x : a 19 x 1 column consisting of the terms, an and bn. Where x^ 

through x are the 10 an's and x through x are the 9 

& 10 11 & 19 

bn # S. 

To solve for the x^'s, the coefficients of the Fourier series, the 
following matrix operation needs to be performed: 

-l 



[\] - [ F ]' 1 [D] 



(4.6.8) . 



The individual elements in F can be determined by performing the 

au 

integration from 0 to it of the expression for - 5 — as presented in 

OX 

i 

section 4.4.1. Similarly , the individual elements in D can be 

aw 

determined by performing the integration of the expression for — as 

o x 

i 

presented in section 4.4.2. To perform the above integrations, as well 
as inverting the F matrix and solving for the x 's, a program was 
written. This "Toroid Program" is listed in Appendix B. 
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The "Toroid Program" was run to compare the results with those 
obtained using B0S0R4. As can be seen from Figures 4.4 and 4.5, the 
results obtained for the displacements are comparable to those obtained 
from B0S0R4 and are of the same general shape. Additionally, when the 
internal energy of the results are compared, the "Toroid Program" is 
only 4.3% greater than that of the B0S0R4 solution. This difference can 
be accounted for by the selection of the functions chosen to represent u 
and v. The functions chosen for u and v only model the displacements 
and do not model the change in displacements or the curvature of the 
toroid . 



4.7 Forces and Moments: 

The next step was to look at the predictions for the forces, N and 

V 

N, , and the moments, M. and M, , and compare them with the results from 
<p u <p 

B0S0R4 . 



The forces and moments predicted by the "Toroid Program" did not 
corollate well with the results obtained from B0S0R4. As can be seen 
from Figures 4.6 and 4.7, the average of the predicted forces is 
approximately the same as the forces obtained from B0S0R4. However, as 
can be seen from Figures 4.8 through 4.11, the predicted moments were 
several orders of magnitude less than those obtained from B0S0R4. 
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NORMAL DISPLACEMENT 



Figure 4.4 

Normal Displacement (u) 
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TAKCff NTIAX DISPLACEMENT 



Figure 4.5 

Tangential Displacement (v) 
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MERIDIAN RESULTANT Ot>«> 



Figure A. 6 



N vs 6 

V 
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CIRCUMF. RESULTANT («»■) 



Figure 4.7 



n . vs e 

4 > 
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MERIDIAN MOMENT (lt>-ln) 



Figure 4.8 



vs e (B0S0R4) 



0.088 
0.08 
0.018 
0.018 
0.014 
0.018 
0.01 
0.008 
0.008 
0.004 
0.008 
0 

- 0.002 
-0.004 
-0.008 
-0.008 
- 0.01 

0 20 40 80 80 100 120 140 180 180 

(DBGRBBS) 



1/4 BUCKLDTG LOAD 




43 



MVRIDIAN MOMENT Ob 



Figure 4.9 

M. vs 8 (Toroid Program) 
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CIRCUMP. MOMENT («» 



Figure 4.10 
M vs 6 (B0S0R4) 
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CmCUWF. MOMENT Ob 
(Tlm«i lOB— 3) 



Figure 4.11 



M vs 9 (Toroid Program) 



1/4- BUCELDTG LOAD 




(DEGRBB3) 



46 



CHAPTER 5 



CONCLUSIONS 

As was shown in Chapter 4, using the assumed functions for the 
displacements in the energy method produced results for the 

displacements that were consistent with those obtained from B0S0R4 . 
However, it was also shown that the "Toroid Program” was very poor at 
predicting the forces and moments of the toroid. 

The differences can be explained by looking at what the assumed 
functions used for the displacements represents. The assumed functions 
for u and v are fairly good models of the displacements but were not 
modeled to meet the slope or curvature of the displacements except at 
the end points, where 9 - 0 and 9 - n. This resulted in good 
predictions for u and v but inaccurate predictions for u' , u' ' and v' , 
all of which are contained in the expressions for the forces and 
moments . 

Some thought has been given to what other assumed functions might 
be used to model the displacements that would more accurately predict 
the forces and moments. It is believed that a complete Fourier series 
representation, instead of just an odd or even series, might better 
predict the slope and curvature of the displacements and therefore the 
forces and moments. To ensure that the complete Fourier series meets 
the boundary conditions, Lagrange multipliers would be required to 
describe the constraint energy (See section 4.2.3). This is an area 
that is left for further investigation. If a complete Fourier Series 
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representation does not adequately predict the forces and moments the 
next step would be to write some kind of finite element or finite 
difference program to solve this problem. 

Since the program written does model the displacements and produces 
results that are consistent with those obtained from B0S0R4 , the "Toroid 
Program” can be used as a preliminary design tool without having to 
resort to a very complex and expensive computer code such as B0S0R4. 
Unfortunately the "Toroid Program” is restricted to analyzing only 
toroids under hydrostatic loading. Complexities in describing the 
geometry of a general shell of revolution will quickly lead the designer 
to the more complex computer codes like B0S0R4. 

It is also believed that the program written could be modified to 
incorporate nonlinear terms and then be used to predict buckling loads. 
This is again an area that will be left for further investigation. 



48 



REFERENCES 



Almroth, B. 0., Sobel, L. H. , and Hunter, A. R. , 1969, "An Experimental 
Investigation of the Buckling of Toroidal Shells", A1AA Journal, Vol 7, 
No. 11, pp. 2185-2186. 

Bowen, J. D. , 1987, "The Elastic Analysis of a Thin Toroidal Shell", 
Naval Engineer Thesis, Massachusetts Institute of Technology, Cambridge, 
MA. 



Bushnell, D. , 1977, "B0S0R4-Program for Stress, Buckling and Vibration 

of Complex Shells of Revolution", STRUCTURAL MECHANICS SOFTWARE SERIES . 
University Press of Virginia, Charlottesville, VA, pp. 11-143. 

Bushnell, D., 1984, "Computerized Analysis of Shells - Governing 

Equations”, Computers & Structures , Vol. 18, No. 3, pp . 471-536. 

Bushnell, D., 1985, COMPUTERIZED BUCKLING ANALYSIS OF SHELLS, Martinus 

Nijhoff, Boston, MA, pp. 41-45. 

Dyck, V. A., Lawson, J. D., Smith, J. A., 1984, FORTRAN /7 7 AN 
INTRODUCTION TO STRUCTURED PROBLEM SOLVING . Reston Publishing Company, 
Inc., Reston, VA, pp. 410-417. 

James, M. L. , Smith, G. M. and Wolford, J. C., 1977, APPLIED NUMERICAL 
METHODS FOR DIGITAL COMPUTATION with FORTRAN and CSMP, Harper & Row, New 
York, NY, pp. 302-312. 

Reissner, E. , 1963, "On Stresses and Deformations in Toroidal Shells of 
Circular Cross Section Which are Acted Upon by Uniform Normal Pressure”, 
Quarterly of Applied Mathematics , Vol. 21, No. 3, pp . 177 - 188. 

Senjanovic' I., 1972, THEORY OF SHELLS OF REVOLUTION . Brodarski 
Institut, Zagreb, Yugoslavia, pp. 174 - 180. 

Shames, I. J. and Dym, C. L. , 1985, ENERGY AND FINITE ELEMENT METHODS IN 
STRUCTURAL MECHANICS . McGraw-Hill, New York, NY, pp. Ill - 128. 

Sobel, L. H. and Flugge , W. , 1967, "Stability of Toroidal Shells under 
Uniform External Pressure", AIAA Journal, Vol 5, No. 3, pp . 425-431. 

Timoshenko, S. and Woinowsky-Krieger , S., 1959, THEORY OF PLATES AND 

SHELLS . McGraw-Hill, New York, NY, pp. 429 - 447. 



49 



APPENDIX A. 



CURVE FITTING PROGRAMS 



The programs listed were used in conjunction with the out put from 
B0S0R4 to determine the coefficients of the assumed Fourier Series. 
Both programs are written in Fortran 77. 

Program 1 determines ten coefficients for an odd Fourier Series 
that is symmetric about 8 equal to n , to fit the data from B0S0R4 for 
the normal displacement of the shell meridian. 

Program 2 determines nine coefficients for an even Fourier Series 
that is anti -symmetric about 8 equal to ir, to fit the data from B0S0R4 
for the tangential displacement of the shell meridian. 



PROGRAM 1 LISTING: 

C CURVE FITTING PROGRAM 
C LEAST SQUARES CURVE FITTING 
C234567 

INTEGER I , J ,M,N 

REAL X(100) , Y( 100) , F(100 ,10) , FT( 10 , 100) , A(10 , 11) , B( 10) 
REAL FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 , F10 , C( 10) 

EXTERNAL FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 , F10 
N - 100 
M - 10 

C READ X-Y VALUES OF DATA POINTS 
DO 10 I-l.N 

READ (UNIT - 2 , FMT - *) X(I) 

READ (UNIT = 3, FMT - *) Y(I) 
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10 CONTINUE 
C GENERATE THE F MATRIX 



DO 20 I-l.N 

F(I , 1) - FI (X(I) ) 

F(I , 2) - F2(X(I) ) 

F(I , 3) - F3(X(I) ) 

F(I,4) - F4(X(I) ) 

F(I , 5) - F5 (X( I ) ) 

F(1 ,6) - F6(X(I) ) 

F(I , 7) - F7 (X(I) ) 

F( 1 , 8) - F8(X(I) ) 

F(I,9) - F9(X(I) ) 

F(I , 10) - F10(X(I) ) 

20 CONTINUE 

C GENERATE THE TRANSPOSE OF THE F MATRIX 
DO 30 1-1 , N 
DO 30 J-l.M 

FT ( J , I ) - F(I,J) 

30 CONTINUE 

C DETERMINE COEFFICIENT MATRIX A OF SIMULTANEOUS 
C EQUATION SYSTEM 

CALL MATMPY (FT,F,A,M,N,M) 

C DETERMINE COLUMN OF CONSTANTS FOR SIMULTANEOUS 
C EQUATION SYSTEM 

CALL MATMPY (FT,Y,B,M,N,1) 

DO 40 I— 1 ,M 
A(I,M+1) - B(I ) 

40 CONTINUE 

C DETERMINE A(n) VALUES BY SOLVING SIMULTANEOUS EQUATIONS 
C USING CHOLE3KY METHOD 
MP1 - M + 1 
CALL CHLSKY(A,M,MP1 ,C) 

C WRITE OUT THE A(n) VALUES 



51 



WRITE (UNIT - 4 , FMT - *) (I ,C(I) , 1-1 ,M) 

END 

C 

C DETERMINES MATRIX C AS PRODUCT OF A AND B MATRICES 
SUBROUTINE MATMPY(A , B , C ,M,N, L) 

REAL A(M,N) ,B(N,L) ,C(M,L) 

DO 10 I-l.M 
DO 10 J-l.L 
C(I,J) - 0. 

DO 10 K-l.N 

C(I,J) - C(I,J)+A(I I K)*B(K,J) 

10 CONTINUE 
END 

C DEFINE THE FUNCTIONS 
C 

C DEFINE THE A(0) FUNCTION 
REAL FUNCTION F1(X) 

REAL X 
FI - 1.0 
END 

C DEFINE THE A(l) FUNCTION 
REAL FUNCTION F2(X) 

REAL X 
F2 - COS(X) 

END 

C DEFINE THE A(2) FUNCTION 
REAL FUNCTION F3(X) 

REAL X 

F3 - COS (2 . *X) 

END 

C DEFINE THE A(3) FUNCTION 
REAL FUNCTION F4(X) 

REAL X 



52 



F4 - COS ( 3 . *X) 

END 

C DEFINE THE A(4) FUNCTION 
REAL FUNCTION F5(X) 
REAL X 

F5 - COS (4 . *X) 

END 

C DEFINE THE A(5) FUNCTION 
REAL FUNCTION F6(X) 
REAL X 

F6 - COS ( 5 . *X) 

END 

C DEFINE THE A(6) FUNCTION 
REAL FUNCTION F7(X) 
REAL X 

F7 - COS ( 6 . *X) 

END 

C DEFINE THE A(7) FUNCTION 
REAL FUNCTION F8(X) 
REAL X 

F8 - COS ( 7 . *X) 

END 

C DEFINE THE A(8) FUNCTION 
REAL FUNCTION F9(X) 
REAL X 

F9 - COS (8 . *X) 

END 

C DEFINE THE A(9) FUNCTION 
REAL FUNCTION F10(X) 
REAL X 

F10 - COS (9 . *X) 

END 

C 
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SUBROUTINE CHLSKY(A,N,M,X) 

REAL A( 10,11) ,X(10) 

C CALCULATE FIRST ROW OF UPPER UNIT TRIANGULAR MATRIX 
DO 10 J-2 ,M 
A( 1 , J) - A(1 , J)/A(l , 1) 

10 CONTINUE 

C CALCULATE OTHER ELEMENTS OF U AND L MATRICES 
DO 60 1-2, N 
J - I 

DO 30 II-J.N 
SUM - 0. 

JM1 - J-l 

DO 20 K-l.JMl 

SUM - SUM+A(II,K)*A(K,J) 

20 CONTINUE 

A( II , J ) - A(II , J) -SUM 
30 CONTINUE 
IP1 - 1+1 
DO 50 JJ-IP1 , M 
SUM - 0. 

IM1 - 1-1 

DO 40 K-l.IMl 

SUM - SUM+A(I ,K)*A(K,JJ) 

40 CONTINUE 

A(I,JJ) - (A(I , JJ) -SUM) /A (I , I) 

.50 CONTINUE 
60 CONTINUE 

C SOLVE FOR X(I) BY BACK SUBSTITUTION 
X(N) - A(N , N+l) 

L - N-l 
DO 80 NN-l.L 
SUM - 0. 

I - N-NN 
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70 



IP1 - 1+1 
DO 70 J-IPl.N 
SUM - SUM+A (I,J)*X(J) 
CONTINUE 

X(I) - A(I ,M) -SUM 
80 CONTINUE 
END 
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PROGRAM 2 LISTING: 



C CURVE FITTING 

C LEAST SQUARES CURVE FITTING 

C234567 

INTEGER I,J,M,N 

REAL X(IOO) , Y(IOO) , F(100 ,9) , FT(9 , 100) , A(9 , 10) , B(9) 
REAL FI , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 , C(9) 

EXTERNAL Fl , F2 , F3 , F4 , F5 , F6 , F7 , F8 , F9 
N - 100 
M - 9 

C READ X-Y VALUES OF DATA POINTS 
DO 10 I-l.N 

READ (UNIT - 2,FMT - *) X(I) 

READ (UNIT - 5 , FMT - *) Y(I) 

10 CONTINUE 
C GENERATE THE F MATRIX 
DO 20 I-l.N 

F( I , 1) - Fl (X( I ) ) 

F(1 , 2) - F2 (X(I) ) 

F(1 , 3) - F3(X(I)) 

F(I ,4) - F4(X(I) ) 

F( I , 5) - F5(X(I) ) 

F(1 , 6) - F6(X(I) ) 

F(1 , 7) - F7 (X(I) ) 

F(I,8) - F8 (X(I) ) 

F(I , 9) - F9(X(I) ) 

20 CONTINUE 

C GENERATE THE TRANSPOSE OF THE F MATRIX 
DO 30 I-l.N 
DO 30 J-l.M 

FT(J.I) - F( I , J ) 

30 CONTINUE 
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C DETERMINE COEFFICIENT MATRIX A OF SIMULTANEOUS 
C EQUATION SYSTEM 

CALL MATMPY (FT,F,A,M,N,M) 

C DETERMINE COLUMN OF CONSTANTS FOR SIMULTANEOUS 
C EQUATION SYSTEM 

CALL MATMPY (FT,Y,B,M,N,1) 

DO 40 I-1,M 
A(I ,M+1) - B(I) 

40 CONTINUE 

C DETERMINE B(n) VALUES BY SOLVING SIMULTANEOUS EQUATIONS 
C USING CHOLESKY METHOD 
MP1 - M + 1 
CALL CHLSKY(A,M,MP1 , C) 

C WRITE OUT THE B(n) VALUES 

WRITE (UNIT - 4 , FMT - *) (I , C( I) , 1-1 ,M) 

END 

C 

C DETERMINES MATRIX C AS PRODUCT OF A AND B MATRICES 
SUBROUTINE MATMPY (A,B,C,M,N,L) 

REAL A(M,N) ,B(N,L) ,C(M,L) 

DO 10 I— 1 ,M 
DO 10 J-l , L 
C(I,J) - 0. 

DO 10 K-l.N 

C(I,J) - C(I,J)+A(I,K)*B(K,J) 

10 CONTINUE 
END 
C 

C DEFINE THE FUNCTIONS 
C 

C DEFINE THE B(l) FUNCTION 
REAL FUNCTION F1(X) 

REAL X 
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FI - SIN (X) 

END 

C DEFINE THE B(2) FUNCTION 
REAL FUNCTION F2(X) 
REAL X 

F2 - SIN(2.*X) 

END 

C DEFINE THE B ( 3 ) FUNCTION 
REAL FUNCTION F3(X) 
REAL X 

F3 - SIN(3.*X) 

END 

C DEFINE THE B(4) FUNCTION 
REAL FUNCTION F4(X) 
REAL X 

F4 - SIN(4 . *X) 

END 

C DEFINE THE B(5) FUNCTION 
REAL FUNCTION F5(X) 
REAL X 

F5 - SIN(5.*X) 

END 

C DEFINE THE B(6) FUNCTION 
REAL FUNCTION F6(X) 
REAL X 

F6 - SIN(6.*X) 

END 

C DEFINE THE B(7) FUNCTION 
REAL FUNCTION F7(X) 
REAL X 

F7 - SIN(7.*X) 

END 

C DEFINE THE B(8) FUNCTION 
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REAL FUNCTION F8(X) 

REAL X 

F8 - SIN (8 . *X) 

END 

C DEFINE THE B ( 9 ) FUNCTION 
REAL FUNCTION F9(X) 

REAL X 

F9 - SIN(9.*X) 

END 

C 

SUBROUTINE CHLSKY(A,N,M,X) 

REAL A(9,10) ,X(9) 

C CALCULATE FIRST ROW OF UPPER UNIT TRIANGULAR MATRIX 
DO 10 J-2 ,M 
A(1,J) - A(1 , J)/A(l , 1) 

10 CONTINUE 

C CALCULATE OTHER ELEMENTS OF U AND L MATRICES 
DO 60 1-2, N 
J - I 

DO 30 II-J.N 
SUM - 0. 

JMl - J-l 

DO 20 K-l.JMl 

SUM - SUM+A ( I I , K ) *A ( K , J ) 

20 CONTINUE 

A(II , J) - A(II , J) -SUM 
30 CONTINUE 
IP1 - 1+1 
DO 50 JJ-IPl.M 
SUM - 0. 

IM1 - 1-1 

DO 40 K-l.IMl 

SUM - SUM+A(I,K)*A(K,JJ) 
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40 CONTINUE 

A(I,JJ) - (A(I , JJ) -SUM) /A (I , I) 

50 CONTINUE 
60 CONTINUE 

C SOLVE FOR X(I) BY BACK SUBSTITUTION 
X(N) - A(N , N+I) 

L - N-l 
DO 80 NN-l.L 
SUM - 0. 

I - N-NN 

IP1 - 1+1 

DO 70 J-IPl.N 

SUM - SUM+A ( I , J ) *X ( J ) 

70 CONTINUE 

X(I) - A(I ,M) -SUM 
80 CONTINUE 
END 
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APPENDIX B. 



TOROID PROGRAM 



This program is set up to determine the ten coefficients of the 
even Fourier Series for the normal displacement, (u) , of the toroid and 
the nine coefficients of the odd Fourier Series for the tangential 
displacement, (v) , of the toroid. 

To determine the shape and magnitude of the normal and tangential 
displacement of the toroid the following equations are used: 

9 

u - £ a cos(n0) 

n - 0 
9 

v - I b sin(n0) 

n 

n - 1 

The program was written in Fortran 77. The subroutines used were 
adopted from those presented by Dyck, Lawson and Smith (1984). 

PROGRAM LISTING: 

C TOROID PROGRAM 
INTEGER I ,N 

REAL F( 19 , 19) ,D(19) ,C(19) 

REAL E , V , RB , RL , H , P , Al , A2 , A3 , A 
LOGICAL ERROR 

EXTERNAL FA, FAl , FB , FBI , SIMP 

PRINT *,'This Program will calculate the coefficients' 

PRINT *,'of the Fourier Series for the Normal and’ 
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PRINT Tangential displacements of a thin, circular' 
PRINT *, 'toroidal shell when loaded by uniform' 

PRINT *, 'external hydrostatic pressure.' 

PRINT * , ' ' 

PRINT *, 'Enter Youngs Modulus, E(psi) -' 

READ * , E 

PRINT *, 'Enter Poissons Ratio, V -' 



READ *,V 

PRINT *, 'Enter Toroid radius, R(in) -' 

READ * ,RB 

PRINT *, 'Enter Circular Section radius, r(in) -' 
READ * ,RL 

PRINT *, 'Enter Shell thickness, h(in) -' 

READ * , H 

PRINT *, 'Enter hydrostatic pressure, P(psi) -' 
READ *,P 

PI - 4 . 0*ATAN (1.0) 

N - 19 
A - RB/RL 

Al-2 . 0*PI*E*H/(1 . -V**2) 



A2--2.0*PI*P*RL**2 



A3- ( - PI*E*(H**3 ) ) / ( 6 . 0* (RL**2 ) * ( 1 . 0 - V**2 ) ) 

C GENERATE THE D COLUMN 
D(1)-A2*PI*A 
D(2)-A2*PI/2.0 
DO 30 1-3, N 
D(I )-0 . 0 
30 CONTINUE 
C GENERATE THE F MATRIX 
DO 40 I-l.N 

PRINT *,' INTEGRATING MATRIX COLUMN', I 

F(1 , I) - Al*SIMP(FA, PI , V , A, I , 0 . 0)+A3*SIMP (FAl , PI, V, A, I, 0.0) 
F(2 , I) - A1*SIMP(FA, PI , V , A, I , 1 . 0)+A3*SIMP(FAl , PI, V, A, I, 1.0) 
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F( 3 , I) - Al*SIMP(FA, PI , V , A , I , 2 . 0)+A3*SIMP(FAl , PI, V, A, I, 2.0) 
F(4, I) - Al*SIMP(FA, PI , V , A, I , 3 . 0)+A3*SIMP(FAl ,PI,V,A,I,3.0) 
F(5 , I) - Al*SIMP(FA, PI , V , A, I , 4 . 0)+A3*SIMP(FAl ,PI,V,A,I,4.0) 
F(6 , I) - A1*SIMP(FA ,PI,V,A,I,5. 0)+A3*SIMP(FAl , PI, V, A, I, 5.0) 
F(7 , I) - A1*SIMP ( FA, PI, V, A, I, 6. 0)+A3*SIMP(FAl , PI, V, A, I, 6.0) 
F(8 , I) - A1*SIMP(FA, PI , V, A, I , 7 . 0)+A3*SIMP(FAl , PI, V, A, I, 7.0) 
F(9 , I) - Al*SIMP(FA t PI,V,A,I,8.0)+A3*SIMP(FAl,PI,V,A,I,8.0) 
F(10 , I) - A1*SIMP(FA, PI , V, A , I , 9 . 0)+A3*SIMP(FAl , PI, V, A, I, 9.0) 
F(ll , I) - A1*SIMP(FB,PI,V,A,I,1.0)+A3*SIMP(FB1,PI,V,A,I,1.0) 
F( 12 , I) - Al*SIMP(FB ,PI,V,A,I,2. 0)+A3*SIMP(FBl ,PI,V,A t I,2.0) 
F(13 , I) - Al*SIMP(FB , PI , V , A, I , 3 . 0)+A3*SIMP(FBl ,PI,V,A,I,3.0) 
F(14 , I) - Al*SIMP(FB , PI , V , A, I , 4 . 0)+A3*SIMP(FBl , PI, V, A, I, 4.0) 
F(15 , I) - Al*SIMP(FB , PI , V, A, I , 5 . 0)+A3*SIMP(FBl ,PI,V,A,I,5.0) 
F( 16 , I) - A1*SIMP(FB,PI 1 V,A,I,6.0)+A3*SIMP(FB1 ,PI i V,A,I, 6.0) 
F( 17 , 1) - Al*SIMP(FB , PI , V , A, 1 , 7 . 0)+A3*SIMP(FBl , PI, V, A, I, 7.0) 
F( 18 , I) - Al*SIMP(FB , PI , V , A, I , 8 . 0)+A3*SIMP( FBI , PI, V, A, 1,8.0) 
F( 19 , I ) - Al*SIMP (FB , PI , V , A , I , 9 . 0)+A3*SIMP(FBl ,PI,V,A, 1,9.0) 
40 CONTINUE 

: CONVERT THE LINEAR SYSTEM TO A TRIANGULAR SYSTEM 

CALL GAUSS (F,D,N, ERROR) 

IF (ERROR) THEN 

PRINT *, 'MATRIX GENERATED IS SINGULAR,' 

PRINT *, 'SOLUTION IS NOT POSSIBLE.' 

ELSE 

C SOLVE THE TRIANGULAR LINEAR SYSTEM 
CALL BSOLVE(F,D,N,C) 

PRINT * , ' THE FOLLOWING ARE THE 10 COEFFICIENTS FOR THE' 



PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT *, ' 



'NORMAL DISPLACEMENT OF THE TOROID.' 



'USE IN THE SERIES; An*COS (n*theta/r) 



n 



An 



n 



An' 
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PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 
PRINT * 



(I-1,C(I), 1-1,10) 

t f 

'THE FOLLOWING ARE THE 9 COEFFICIENTS FOR THE' 
'TANGENTIAL DISPLACEMENT OF THE TOROID.' 

t t 

'USE IN THE SERIES; Bn*SIN(n*theta/r) ' 



n 



Bn 



n 



+ 



Bn' 



PRINT *, (I-IO.C(I) ,1-11,19) 

END IF 

END 



C 

C THIS SUBPROGRAM FUNCTION INTEGRATES THE MATRIX FUNCTIONS 
C DEFINED BY USING SIMPSON'S RULE AS AN APPROXIMATION 
REAL FUNCTION SIMP(F,PI ,V,A,N,Z) 

REAL PI , F , H , SUMEVN , SUMODD , X , B 

INTEGER I 

EXTERNAL F 

B-0.0 

H-PI/100 . 

SUMEVN— 0 . 0 
SUMODD-F ( A , V , N , Z , H) 

DO 100 1-1,49 
X-2 . *I*H 

SUMEVN-SUMEVN+F ( A , V , N , Z , X) 

SUMODD-SUMODD+F ( A , V , N , Z , X+H) 

100 CONTINUE 

SIMP-(H/3 .)*(F(A,V,N,Z,B)+4. *SUMODD+2 . *SUMEVN+F(A, V ,N , Z , PI ) ) 
RETURN 
END 
C 

C 
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C DEFINE THE MATRIX FUNCTIONS 
C CALCULATE ROWS 1 THROUGH 10 
C FUNCTION FA REPRESENTS MEMBRANE ENERGY 
REAL FUNCTION FA(A,V, I , Z,X) 

REAL Q 

IF (I .LT. 11) THEN 
Q-I-l 

FA— (A+COS (X) )*COS (Q*X)*COS (Z*X)+2 . 0*V*COS (X)*COS (Q*X)* 
+ COS ( Z*X)+(COS (X) **2 )*COS (Q*X)*COS ( Z*X)/(A+COS (X) ) 

ELSE 

Q-I-10 

FA- (A+COS (X) )*Q*COS (Q*X)*COS (Z*X)+V*(COS (X)*Q*COS (Q 
+ *X)*COS(Z*X) -SIN(X)*SIN(Q*X)*COS(Z*X)) -COS(X)*SI 

+ N(X)*SIN(Q*X)*COS (Z*X)/( A+COS (X) ) 

END IF 
END 

C FUNCTION FAl REPRESENTS BENDING ENERGY 
REAL FUNCTION FAl(A,V, I ,Z,X) 

REAL Q 

IF (I .LT. 11) THEN 
Q-I-l 

FAl-(A+COS(X))*((Q*Z)**2)*COS(Q*X)*COS(Z*X) -V*SIN(X 
+ )*(Q*(Z**2)*SIN(Q*X)*COS(Z*X)+(Q**2)*Z*COS (Q*X)* 

+ SIN(Z*X) )+(SIN(X)**2)*Q*Z*SIN(Q*X)*SIN(Z*X)/(A+C 

+ OS(X)) 

ELSE 

Q-I-10 

FAl- (A+COS (X) )*Q*(Z**2)*COS (Q*X)*COS (Z*X) -V*SIN(X)* 

+ (Q*Z*COS (Q*X)*SIN(Z*X)+(Z**2)*SIN (Q*X)*COS (Z*X) ) 

+ +( SIN (X)**2)*Z*SIN(Q*X)*SIN(Z*X)/ (A+COS (X) ) 

END IF 
END 
C 
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C CALCULATE ROWS 11 THROUGH 19 
C FUNCTION FB REPRESENTS MEMBRANE ENERGY 
REAL FUNCTION FB(A, V , I , Z ,X) 

REAL Q 

IF (I .LT. 11) THEN 
Q-I-l 

FB-(A+COS (X) ) *Z*COS (Q*X) *COS (Z*X)+V* (COS (X)*Z*COS(Q 
+ *X)*COS (Z*X) -SIN(X)*COS (Q*X)*SIN(Z*X) )-COS(X)*SI 

+ N(X)*COS(Q*X)*SIN(Z*X)/ (A+COS(X) ) 

ELSE 

Q-I-10 

FB— (A+COS (X) )*Q*Z*COS (Q*X)*COS ( Z*X) - V*SIN(X>* ( Z*SIN 
+ (Q*X)*COS (Z*X)+Q*COS (Q*X)*SIN( Z*X) )+(SIN(X)**2)* 

+ SIN(Q*X)*SIN(Z*X)/(A+COS(X)) 

END IF 
END 

C FUNCTION FBI REPRESENTS BENDING ENERGY 
REAL FUNCTION FB1(A,V,I,Z,X) 

REAL Q 

IF (I .LT. 11) THEN 
Q-I-l 

FBI- (A+COS (X) )*Z*(Q**2)*COS (Q*X)*COS (Z*X) -V*SIN(X)* 
+ (Q*Z*SIN(Q*X)*COS (Z*X)+(Q**2 )*COS (Q*X)*SIN(Z*X) ) 

+ +( S IN (X)**2)*Q*SIN(Q*X)*SIN(Z*X)/ (A+COS (X) ) 

ELSE 

Q-I-10 

FBI— (A+COS (X) ) *Z*Q*COS (Q*X) *COS ( Z*X) - V*SIN(X)*(Z*SI 
+ N(Q*X)*COS (Z*X)+Q*COS (Q*X)*SIN(Z*X) )+(SIN(X)**2) 

+ *SIN(Q*X)*SIN(Z*X)/( A+COS (X) ) 

END IF 
END 
C 

C GAUSSIAN ELIMINATION MODULE 
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SUBROUTINE GAUSS (A, B ,M, ERROR) 
INTEGER M, I , J ,K, INDEX 
REAL A(M,M) ,B(M) .PIVOT, TEMP, RATIO 
LOGICAL ERROR 
ERROR- . FALSE . 

DO 100 I-l.M 
INDEX-I 

PIVOT— ABS ( A( I , I ) ) 

DO 200 J-I+l.M 

IF (ABS (A( J , I) ) .GT. PIVOT) THEN 
PIVOT-ABS ( A( J , I ) ) 

INDEX-J 
ENDIF 
200 CONTINUE 

IF (INDEX .GT. I) THEN 
DO 400 K-I.M 
TEMP-A(I,K) 

A(I .K)-A(INDEX.K) 

A ( INDEX, K) -TEMP 
400 CONTINUE 
TEMP-B(I) 

B(I)-B (INDEX) 

B( INDEX) -TEMP 
ENDIF 

IF (PIVOT .EQ. 0.0) THEN 
ERROR - .TRUE. 

ELSE 

DO 300 J-I+l.M 

RATIO-A(J ,I)/A(I,I) 

DO 500 K-I+l.M 

A(J ,K)-A(J ,K) -A(I ,K)*RATIO 
500 CONTINUE 

B( J)-B(J) -B(I)*RATIO 
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300 CONTINUE 
ENDIF 

100 CONTINUE 
END 
C 

C SUBROUTINE TO SOLVE TRIANGULAR LINEAR SYSTEM 
SUBROUTINE BSOLVE(A , B ,M, Z) 

INTEGER M , I , J 
REAL A(M,M) ,B(M) ,Z(M) , SUM 
DO 200 I— M, 1,-1 
SUM-B(I) 

DO 100 J-I+l.M 

SUM-SUM-A(I , J)*Z(J) 

100 CONTINUE 

Z ( I ) -SUM/A (1,1) 

200 CONTINUE 
END 
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